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Abstract 

A numerical approach to ground-state dynamical correlation functions from 
Density Matrix Renormalization Group (DMRG) is developed. Using sum 
rules, moments of a dynamic correlation function can be calculated with 
DMRG, and with the moments the dynamic correlation function can be ob- 
tained by the maximum entropy method. We apply this method to one- 
dimensional spinless fermion system, which can be converted to the spin 1/2 
Heisenberg model in a special case. The dynamical density-density correlation 
function is obtained. 

PACS Numbers: 71.10.+X, 75.10.Jm, 75.40.Gb 
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Dynamical correlation functions of a model are of special interest, because they can pro- 
vide a comprehensive comparison to experimental measurements. Unfortunately they are 
very difficult to calculate analytically or numerically for strongly correlated systems. Even 
for one-dimensional systems, the dynamical correlation functions are hard to obtain. For 
example the 5 = 1/2 Heisenberg model, although its exact solution from Bathe Ansatz 
has been known for a long time, its ground-state dynamical correlation functions have not 
yet been obtained. Until now there are only a few general ways to obtain dynamical cor- 
relation functions. Analytically, only the asymptotic behavior of correlation functions for 
one-dimensional models in the quantum critical regime are able to obtain by bosonization 
or conformal field theory Jl| . Numerically, one way to calculate dynamical correlation func- 
tions is the analytic continuation of quantum Monte Carlo simulations with the maximum 
entropy method @J3[]. But this method will encounter essential difficulties if we are interested 
in zero-temperature properties. Another numerical method to calculate the ground-state dy- 
namical properties [|J is based on the Lanczos method. The limitation of this method is 
that it cannot be applied to large systems. 

The Density Matrix Renormalization Group (DMRG) method proposed by White || 
is a powerful method to study the ground state of one- dimensional interacting systems. 
With this method the ground state energy, a few excitation energies, and static correlation 
functions can be calculated for a large system. However, it was not clear if one can obtain 
dynamical properties from this method. 

In this paper, we describe a numerical method for calculating ground-state dynamical 
correlation functions in a systematic way, which is a combination of DMRG and maximum 
entropy methods (MEM) ||. In general the moments of a dynamical correlation function 
can be expressed as static correlation functions, which can be calculated by DMRG method. 
With these moments we can obtain the dynamic correlation function with MEM. We apply 
this method to the one-dimensional spinless fermion system with nearest neighbor interac- 
tion. This model is equivalent to the spin-1/2 XX Z chain. We have considered two special 
cases of this model, corresponding to the XY model and the Heisenberg model. The dy- 
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namical density-density correlation (namely the structure function S(q, u) in spin chain) is 
obtained. For non-interacting case (the XY model) we compare our result with the exact 
result, and obtain a very good agreement. 

The one-dimensional spinless fermion model we consider has the following Hamiltonian: 

H = -t ^(cjci+i+h.c.) + V Ei ^n i+ i, (1) 

i 

where cjp are annihilation (creation) operators for a fermion at site i, and rti = c\ci — |. 
The Hamiltonian written in such form ensures the ground state is at half filling. This model 
may be mapped to the XXZ model by the Jordan- Wigner transformation. Under this 
transformation Sf = Hi, J x = J y = 2t, and J z = V. At V = this model is equivalent to 
the XY model, while at V = 2t it is equivalent to the Heisenberg model. We only consider 
these two cases in this paper, and the results for other V will be presented elsewhere [0. 

The first step of our method is to use sum rules to express the moments of a dynamical 
correlation function by some static correlation functions. The sum rules for the spin model 
have been derived IHl. We use the similar definition of the correlation functions as in Ref. 



Xc(q,t) = ^({n(q,t),n(-q,0)}) - (n(q,t))(n(-q,0)), 

X"(q,t) = ^({n(q,t),n(-q,0)]), (2) 

where n(q) = N~ 1 / 2 J2 n i eiql i the curly bracket is an anticommutator, and 
(n(q)) =Tr{n{q)e~" ). The fluctuation-dissipation theorem gives the relation: Xc{q, = 
coth(u /2k B T)x" (q, oS). The structure function or dynamic form factor S(q,u) is defined 
as S(q,u) = x"(g, u;)/(l — e~ w / fcsT ). Due to the parity and time reversal symmetry in 
our model, x"(<2S ^ an d Xc(q,u) have following properties: x"(<?> ~ u ) = ~x"( ( l, UJ ) an d 
Xc{q, = Xili^)- At zero temperature S(q,ui) = Xc(q,u) = x"{Q.i u ) f° r u > 0, therefore 
the sum rules given in Ref. || can be written as 

r°° duj X "{q,uj) 1 

mi{q)= = -x{q,uJ = 0) 

Jo it to 2 



m 2 (q) 
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duj x"(q,v) 



Xc(q,t = 0) 



m 3 (q) 



f 

Jo 



oo 




i<[[^n(g)],n(-g)]> 



2(cjcj + i)(l - cos(g)) 



(3) 



where x(<L w = 0) is the static susceptibility. These sum rules can be easily generalized to 
higher moments: 



Apart from the first moments which is given by the static susceptibility, all the other mo- 
ments can be expressed as equal-time correlation functions. Theoretically if all the moments 
are known, one can obtain the x"(q,t) an d thus x"^^)- 111 rea l calculations, it is tedious 
to calculate the commutators for higher moments, and there are more and more new equal- 
time correlation functions appear in the expression of higher moments. However it is still 
reasonable to obtain the expression for the first several moments using a symbolic manipu- 
lator, such as Mathematica, to calculate the commutators. In this work we have calculated 
the expressions for the first five moments. Details of the expressions for the fourth and fifth 
moments will be given elsewhere 

The second step is to obtain the moments by calculating those static correlation functions 
with DMRG. The infinite lattice method (see Ref. || for details) is used in our calculations 
for open ended chains, t — 1 is chosen, and states kept at each iteration varies from 52 to 64. 
We calculate the equal-time correlations, for example (niUj), by taking j in the middle of the 
system. For a system which has parity and translational symmetries, (rijUi) only depends 
on \i — j\. Therefore (n g n_ q ) = J2i( n j n j+i) etql is independent of j. Since the calculations 
are done with open boundary condition, (n J n J+ ;) depends on the position j. The boundary 
effect is larger when j or j + I is closer to boundary, therefore we choose j at the center 





-±([[H,...,[H,n(q)]...},n(-q)}) I even 
l{{[H,...,[H,n(q)]...],n(-q)}) I odd 
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of the system. Also due to the open boundary, the correlation (njUj+i) has an even-odd 
oscillation in j. We take the mean value of (rijUj+i) at j even and odd, which is close to the 
value with period boundary condition for a system having the same size. When the system 
size goes to infinity, the boundary effect can be neglected. We calculate the moments for 
system sizes varying from 100 sites to 200 sites, and obtain their values for infinite system 
by extrapolating the data. 

The next step is to use MEM to obtain the dynamical correlation functions. MEM 
has become standard way to extract maximum information from incomplete data || . This 
method has been applied to the analytic continuation of the Quantum Monte Carlo data 
0, and in this paper we apply a similar method to extract the dynamic susceptibility 
~x"{q,oS) from the finite number of moments mi with the corresponding errors 0\. Defining 
f(u>) = x"(u)/uj, which is a positive definite, as the distribution function, and the entropy 
or the information function S = £) w f{ u ) ~ fi 00 ) log/C^)- By maximizing the entropy under 
the constrains mi — J °° —u) l ~ x f{uS) = 0, f{uj) has the following form 

m = e -sr= 1 (w- i ) ) (4) 

where n is the number of moments and A; are the Lagrange multipliers. At this point one may 
try to find A; by requiring the f(ui) to satisfy the constrains without considering the error bars 
of the moments. However, in general, the error bars cannot be neglected. The kernel of the 
transformation is singular, so small error bars in moments may produce large errors in f(uj). 
By maximizing the posterior probability e aS ~ L where L = Y,i( m i ~ Io° ujl ~ 1 f{ UJ )) 2 / (T ij one 
can find the most probable f(u), which gives us the moments within the range of error bars. 

Let us first discuss the extrapolation and the error bar of our DMRG results. There are 
two major contributions to the error: that from finite size effects and that from basis set 
transaction in the DMRG calculations. The error bar of DMRG calculation for any finite 
size is obtained by varying the number of states kept at each iteration, whereas the finite 
size error is obtained by varying the system size. The asymptotic behavior of correlation 
functions is known for this model , which decay as a power of system size. For a system 
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with a gap, the extropolation should be done as an exponential function of system size. 
In Fig. [I], we plot the second and third moments at q = n/2 for V = 2t as a function of 
1/N, where N is the number of sites of the system. The error from basis set transaction 
produces the error in the extrapolating values. We use this resultant error to estimate the 
error bar of the moments. Extrapolating to 1/N — > gives m 2 = 0.1700, and the error bar 
is estimated as 10~ 4 . For the third moment we have = 0.59085 and the estimated error 
bar is 2 x 10 -5 . Actually the third moment is known exactly: m% = — |i?(l — cos(g)) with 
the ground state energy per site E = — 2 (In 2 — 1/4). The exact value at 7r/2 is 0.590863. 
We test our method for the non-interacting case. In this case, )("{q,uj) is known exactly: 

6{uj-2t\sm{q)\)d{tt\ S m{q/2)\-u) 
X W ' } [mHm 2 {q/2)-uY 12 ' [ > 

where 9(x) is the step function. The moments can also be calculated analytically. In Table 
|, we compare the moments calculated by DMRG with the exact results. The error bars 
obtained by DMRG provide a very good estimation. Apart from the five moments, there 
are two more pieces of information in this case: the energy boundaries 2t|sin(g)| < uj < 
4t| sin(g/2)| for x"(q, u )- Using the MEM, we obtain x(g, uj) for q = 2n /3. In Fig. 0, we plot 
x{q, uj) obtained by MEM with different number of moments and the exact one from Eq. 
(5). It shows that the x(<?, ^) obtained by MEM converge to the exact one when the number 
of moments is increased, and x(g, u) calculated with five moments is a good approximation 
for the exact result. We have also calculated x(g, uj) for other q, they have similar behavior. 

For the interacting case with V = 2t, which corresponding to the Heisenberg model, the 
elementary excitations are known as 5* = 1/2 objects |10[ (spinons). The dispersion relation 
is e(q) = 7rt|sin(g)| |TTJ, which provides the lower bound of excitation energies for each 
momentum q. The spectral weight is dominated by the continuum of the two-spinon excited 



states |12[], and the energy range for the continuum is tir\ sin(g)| < uj < 2t|7rsin(g/2)|. Since 
the contributions from the excited states of more than two spinons are finite, we only have 
the low energy bound. In Fig. ^ the x"( ( ?! to ') obtained by MEM with different number of 
moments are plotted for q = 2tt/3. One can see the tendency of the curves as the number of 
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moments increase. x"(Qj a; ) tends to diverge at lower-bound. In Fig. f|, x"{q.i u ) * s plotted for 
other momentums. We marked the position of upper-bound for two-spinon excited states. It 
is obvious that the contributions from excited states with more than two spinons are finite, 
although they are small. 

In conclusion, we have developed a numerical method for calculating the ground-state 
dynamical correlation functions in one- dimensional quantum systems based on the Density 
Matrix Renormalization Group Method and the maximum entropy method. We demonstrate 
this method on the dynamical density-density correlation x"{li u ) °f the spinless fermion 
system with nearest neighbor interaction. For the non-interacting case, it corresponds to the 
XY model, and the dynamical density-density correlation function obtained by our method 
shows a very good agreement with the exact result. For the interacting case with V = 2t, it 
corresponds to the Heisenberg model, we obtain the x"(li u )i which was not known before. 
This method is a very general one, which can be applied to any one-dimensional system 
with short range interaction like, e.g. the Hubbard model, the S = 1 Heisenberg model, the 
interacting fermion (or boson) system with randomness. 

We would like to acknowledge useful discussions with J.E. Gubernatis, Shoudan Liang, 
and R.N. Silver. This work was supported by the National Science Foundation grant No. 
DMR-9107563. In addition MJ would like to acknowledge the support of the NSF NYI 
program. 
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FIG. 1. Moments versus the inverse system size 1/N — > for V = 2t and q = it/ 2, where 
M is the number of states kept at each iteration in DMRG calculations. The extropolation is to 
1/N — > 0, and the error is estimated by the different extropolations caused by the errors in slope. 
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FIG. 2. The dynamical structure function x"(<Z) a; ) for F = (the XY model) and q = 2ir/?>. 
We plot the results obtained by MEM with different number of moments. Two solid vertical lines 
are the energy boundaries. 
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FIG. 3. The dynamical structure function x"{Qi u ) f° r V = 2t (the Heisenberg model) and 
q = 27r/3. We plot the results obtained by MEM with different number of moments. The solid 
vertical line is the lower boundary. The arrow marks the position of the upper boundary for the 
two-spinon excited states. 
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FIG. 4. The dynamical structure function x"(c/, u->) for V = 2t (the Heisenberg model) at (a) 
q = it/2, and (b) tt/3 and tt/4. x'^Q, 1 ^) are obtained by MEM with five moments. The solid 
vertical lines are the lower boundaries for each momentum. The arrow marks the position of the 
upper boundary of the two-spinon excited states. 
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TABLES 

TABLE I. The comparison of the moments obtained by DMRG with the exact results for V = 

and q = 2tt/3. In the DMRG calculations 64 states are kept at each iteration, and the error bars 
are estimated by changing the number of states kept and finite size scaling. 

EXACT DMRG ERROR 

m 1 (2vr/3) 0.121013 0.1211 1 x 10~ 3 

m 2 (2^/3) 0.333333 0.33337 5 x 10~ 5 

m 3 (2vr/3) 0.954930 0.954928 5 x 10~ 6 

m 4 (27r/3) 2.826993 2.8273 5 x 10~ 4 

m 5 (2^/3) 8.594367 8.59434 5 x 10~ 5 
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